install.packages("spdep")
library(spdep)
install.packages("spatialreg")
library(spatialreg)
install.packages("tidyr")
library(tidyr)
install.packages("dplyr")
library(dplyr)
install.packages("psych")
library(psych)
install.packages("haven")
library(haven)
install.packages("readxl")
library(readxl)
install.packages("writexl")
library(writexl)

#set WD
#download data file 'medata_1975_2021'

medata_1975_2021 <- read_excel("medata_1975_2021.xlsx")

colnames(medata_1975_2021)
#[1] "ccode"       "year"        "country"     "me_abs"      "gdp_abs"     "pop"         "nato"        "sample"     
#[9] "StateAbb"    "StateNme"    "yearF"       "ccodeF"      "me_abslag"   "gdp_abslag"  "pop_lag"     "lngdp_abs"  
#[17] "lnpop"       "lnme_abs"    "lnme_abslag" "Lru_abs"     "dist_ru"     "lLru_abs"  

#
#Table 1: Descriptive statistics 
#
#lnme_abs = Milsp (logged) 
#lngdp_abs = GDP (logged)
#lnpop = Pop (logged) 
#lLru_abs = Russp_(t-1) (logged) 
#dist_ru = Dist

ds1 = describe(medata_1975_2021[, c("lnme_abs", "lngdp_abs", "lnpop", "lLru_abs", "dist_ru")]) #using data file 'medata_1975_2021'
print(ds1)
#vars   n    mean      sd  median trimmed    mad    min     max   range  skew kurtosis    se
#lnme_abs     1 897    8.84    1.92    8.71    8.84   1.65   4.37   13.68    9.31  0.12     0.08  0.06
#lngdp_abs    2 897   26.61    1.65   26.46   26.60   1.86  23.03   30.64    7.61  0.12    -0.45  0.05
#lnpop        3 897   16.53    1.52   16.22   16.61   1.84  12.79   19.62    6.83 -0.33    -0.19  0.05
#lLru_abs     4 652   10.61    0.43   10.68   10.64   0.47   9.59   11.21    1.63 -0.59    -0.61  0.02
#dist_ru      5 652 2575.34 1667.00 2328.00 2246.46 862.87 603.00 7822.00 7219.00  2.16     4.21 65.28

#hns_pctgdp3 = Host-nation support (HNS) >> using data file 'ReplicationData_merged'

#download data file 'ReplicationData_merged'
ReplicationData_merged <- read_excel("ReplicationData_ merged.xlsx")
ds2 = describe(ReplicationData_merged[, "hns_pctgdp3", drop = TRUE]) #option 'drop' converts col into a vector 
print(ds2)
#   vars   n mean   sd median trimmed  mad min max range skew kurtosis se
#X1    1 118 0.02 0.03      0    0.01 0.01   0 0.1   0.1 1.86     2.32  0

ds = rbind(ds1, ds2)
ds.df = as.data.frame(ds)
list.vars = c("Milsp (logged)","GDP (logged)","Pop (logged)","Russp_(t-1) (logged)","Dist","Host-nation support (HNS)")
ds.df$vars = list.vars
colnames(ds.df)
ds.df = ds.df[,-5:-7]
ds.df = ds.df[,-7:-10]
list.col = c("Variable", "N", "Mean", "SD", "Min", "Max")
colnames(ds.df) = list.col

write_xlsx(ds.df, "Table1_SM.xlsx")

#
#Tables 2 and 3
#

#SAC models (1975-2021) in Tables 2 and 3 (in Suppl materials)
#####################################################################################

#based on data file 'medata_1975_2021'
#medata_1975_2021 <- read_excel("medata_1975_2021.xlsx")

data <- medata_1975_2021
data.ordered <- data[order(data$year, data$ccode),] #to be sure of appropriate ordering
#View(data.ordered)

#download a corresponding matrix (inverse distances)
## matrices differ because of unbalanced panel data  
mematrix1975_2021 <- read_excel("mematrix1975_2021.xlsx")
#View(mematrix1975_2021)
capdistinv.wide <- as.matrix(mematrix1975_2021)

capdistinv.2listw.nors <- mat2listw(capdistinv.wide) #'spdep' package;
                                                     ##non-normalized weights matrix, therefore, 'Warning message'  
#class(capdistinv.2listw)

capdistinv.wide.resc <- capdistinv.wide*100 #re-scale to inverse distances in hundreds of km
capdistinv.2listw.nors.resc <- mat2listw(capdistinv.wide.resc) #non-normalized re-scaled weights matrix

counter=seq(1,47,1) #time (trend) counter
y=seq(1975,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered=merge(data.ordered, counter, by="year")
data.ordered$counter.sq=data.ordered$counter^2
data.ordered$counter.cub=data.ordered$counter^3

data.ordered = data.ordered[order(data.ordered$year, data.ordered$ccode),] #to be sure of appropriate ordering for spatial analysis

#Model 1 (Table 2 in Suppl material: Static spatial autocorrelation (SAC) models) 
##without LDV 
sac.resc <- sacsarlm(lnme_abs ~ lngdp_abs + lnpop + as.factor(ccode)+counter+counter.sq+counter.cub, data = data.ordered, 
                     capdistinv.2listw.nors.resc, Durbin = F) #weights non-normalized; 'spatialreg' package
sac.resc$rho

#Model 1 (Table 3 in Suppl material: Dynamic spatial autocorrelation (SAC) models)
##with LDV (spatio-temporal model)
sacL.resc <- sacsarlm(lnme_abs ~ lnme_abslag  + lngdp_abs + lnpop + as.factor(ccode)+counter+counter.sq+ 
                        counter.cub,data = data.ordered, capdistinv.2listw.nors.resc, Durbin = F)
sacL.resc$rho
#summary(sacL.resc)

#####################################################################################

#SAC models (1991-2021) in Tables 2 and 3 (in Suppl materials)
#####################################################################################

data.ordered1991 = subset(medata_1975_2021, medata_1975_2021$year > 1990)
data.ordered1991 <- data.ordered1991[order(data.ordered1991$year, data.ordered1991$ccode),] #to be sure of appropriate ordering for spatial analysis

#download corresponding matrix (1991-2021)
mematrix1991_2021 <- read_excel("mematrix1991_2021.xlsx")
capdistinv.wide1991 <- as.matrix(mematrix1991_2021)

capdistinv.wide1991.resc <- capdistinv.wide1991*100 #re-scale to inverse distances in hundreds of km
capdistinv.2listw.nors.resc1991 <- mat2listw(capdistinv.wide1991.resc) #non-normalized, hece, 'Warning message'

counter=seq(1,31,1)
y=seq(1991,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered1991=merge(data.ordered1991, counter, by="year")
data.ordered1991$counter.sq=data.ordered1991$counter^2
data.ordered1991$counter.cub=data.ordered1991$counter^3

data.ordered1991 = data.ordered1991[order(data.ordered1991$year, data.ordered1991$ccode),] #to be sure of appropriate ordering

#Model 2 (Table 2 in Suppl materials)
##without LDV
sac.resc1991 <- sacsarlm(lnme_abs ~ lngdp_abs + lnpop + as.factor(ccode) +counter+counter.sq+counter.cub, data = data.ordered1991, capdistinv.2listw.nors.resc1991, Durbin = F) #weights non-normalized
sac.resc1991$rho
#summary(sac.resc1991)

#Model 2 (Table 3 in Suppl materials)
##with LDV
sacL.resc1991 <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + as.factor(ccode)+counter+counter.sq+counter.cub, data = data.ordered1991, capdistinv.2listw.nors.resc1991, Durbin = F) #weights non-normalized
sacL.resc1991$rho
#summary(sacL.resc1991)
#####################################################################################

#SAC models (1999-2021) in Tables 2 and 3 (in Suppl materials)
#####################################################################################

data.ordered1999 = subset(medata_1975_2021, medata_1975_2021$year > 1998)
data.ordered1999 <- data.ordered1999[order(data.ordered1999$year, data.ordered1999$ccode),] #to be sure of appropriate ordering for spatial analysis

#download corresponding matrix (1999-2021)
mematrix1999_2021 <- read_excel("mematrix1999_2021.xlsx")
capdistinv.wide1999 <- as.matrix(mematrix1999_2021)

capdistinv.wide1999.resc <- capdistinv.wide1999*100 #re-scale to inverse distances in hundreds of km
capdistinv.2listw.nors.resc1999 <- mat2listw(capdistinv.wide1999.resc) #non-normalized, hence, 'Warning message'

#counter
counter=seq(1,23,1)
y=seq(1999,2021,1)
counter=as.matrix(cbind(y, counter))
colnames(counter) = c("year","counter")
data.ordered1999=merge(data.ordered1999, counter, by="year")
data.ordered1999$counter.sq=data.ordered1999$counter^2
data.ordered1999$counter.cub=data.ordered1999$counter^3

data.ordered1999 = data.ordered1999[order(data.ordered1999$year, data.ordered1999$ccode),] #to be sure of appropriate ordering

#Model 3 (Table 2 in Suppl materials)
sac.resc1999 <- sacsarlm(lnme_abs ~ lngdp_abs + lnpop + as.factor(ccode)+counter+counter.sq+counter.cub, data = data.ordered1999, capdistinv.2listw.nors.resc1999, Durbin = F)
sac.resc1999$rho
#summary(sac1999.resc)

#Model 3 (Table 3 in Suppl materials)
sacL.resc1999 <- sacsarlm(lnme_abs ~ lnme_abslag + lngdp_abs + lnpop + as.factor(ccode)+counter+
                            counter.sq+counter.cub, data = data.ordered1999, capdistinv.2listw.nors.resc1999, Durbin = F) #weights non-normalized
sacL.resc1999$rho
#summary(sacL.resc1999)

#####################################################################################

###Output
#install.package("texreg")
#library(texreg)

#Models 1-3 in Table 2 in Suppl materials (spatial models)

capture.output(screenreg(list(sac.resc, sac.resc1991, sac.resc1999), stars = c(0.1, 0.05, 0.01), digits=3), file = "Table2_SM.txt")

#Models 1-3 Table 3 in Suppl materials (spatio-temporal models)

capture.output(screenreg(list(sacL.resc, sacL.resc1991, sacL.resc1999), stars = c(0.1, 0.05, 0.01), digits=3), file = "Table3_SM.txt")


####################################################################################
#
# Effects (SR and LRSS in Models 1-3, Tables 2 and 3 in Suppl materials)
#
####################################################################################

#SR for spatial models (Models 1-3 in Table 2 in Suppl materials)
####################################################################################

matrix27x27 <- capdistinv.wide[871:897,871:897] #1-period spatial matrix of 27 NATO member states
matrix27x27.resc <- capdistinv.wide.resc[871:897,871:897] 
dim(matrix27x27.resc)

I <- diag(27) # create identity matrix (N = 27)

#1975
fx <- solve(I - (coef(sac.resc)[1]*matrix27x27.resc)) #1975-2021
#View(fx)
indirect1975 <- mean(rowSums(fx)-diag(fx)) 
indirect1975
#indirect1975 <- mean(colSums(fx)-diag(fx)) 

#1991
fx <- solve(I - (coef(sac.resc1991)[1]*matrix27x27.resc)) #SAC
#View(fx)
indirect1991 <- mean(rowSums(fx)-diag(fx)) 
indirect1991
#indirect1991 <- mean(colSums(fx)-diag(fx)) 

#1999
fx <- solve(I - (coef(sac.resc1999)[1]*matrix27x27.resc)) #1999-2021; (I - rho*W)^(-1)
#View(fx)
indirect1999 <- mean(rowSums(fx)-diag(fx)) 
indirect1999
#indirect1999 <- mean(colSums(fx)-diag(fx)) 

list.indirect = c(indirect1975, indirect1991, indirect1999)
m1 = matrix(NA, 1, 4)
colnames(m1) <- c("", "1975-2021", "1991-2021", "1999-2021")
m1[,1] = "indirect"
m1[1,2:4] = list.indirect
m1.df=as.data.frame(m1)

write_xlsx(m1.df, "Table2_SR spatial effects_SM.xlsx")

####################################################################################


#SR and LRSS effects for spatiotemporal models (Models 1-3 in Table 3 in Suppl materials)
##SR and LRSS effects (1975-2021) for Model 1 (Table 3)
####################################################################################
#
#building L multiplier (for unbalanced data)
#1975-2021
#
#for unbalanced panel data, we have to build L multiplier 

#1975-1981: 13 NATO member states(no data on Canada's gdp until 1996 (incl))
L1 <- diag(13*6) 
TL1 <- rbind((matrix(0,13,dim(L1)[2])),L1) #the firs "layer" of rows represents year 1975 (all 0, hence a matrix of 0)
TL1.t <- cbind(TL1,(matrix(0,dim(TL1)[1], (897 - (dim(TL1)[2])))))

#1982-1996: 13 + Spain
#1982
L2_1 <- matrix(0,14,13) #diag matrix for 1982
for(i in 1:6){
  L2_1[i,i] <- 1
}
for(i in 8:14){
  L2_1[i,(i-1)] <- 1
}
TL2_1 <- cbind((matrix(0,14,(dim(TL1)[2]))),L2_1)
TL2_1.t <- cbind(TL2_1,(matrix(0,dim(TL2_1)[1],(897 - dim(TL2_1)[2]))))

#1983-1996
L2_2 <- diag(14*14) #(14 countries)*(14 years)
TL2_2 <- cbind((matrix(0,dim(L2_2)[1],dim(TL2_1)[2])),L2_2)
TL2_2.t <- cbind(TL2_2, matrix(0,dim(TL2_2)[1],(897 - dim(TL2_2)[2])))

#1997-1998: 15 (with Canada)
L3_1 <- matrix(0,15,14)
L3_1[1,1] <- 1 #US
for(i in 3:15){
  L3_1[i,(i-1)] <- 1
}
TL3_1 <- cbind((matrix(0,15,dim(TL2_2)[2])), L3_1)
TL3_1.t <- cbind(TL3_1, (matrix(0,dim(TL3_1)[1],(897 - dim(TL3_1)[2]))))

#year 1998
L3_2 <- diag(15)
TL3_2 <- cbind((matrix(0,15,dim(TL3_1)[2])), L3_2)
TL3_2.t <- cbind(TL3_2, (matrix(0,dim(TL3_2)[1], (897-dim(TL3_2)[2]))))


#1999-2003
#1999
L4_1 <- matrix(0,18,15) #18 - with 3 additional members, 15 - from the previous year
for(i in 1:10){
  L4_1[i,i] <- 1
}
for(i in 14:18){
  L4_1[i,(i-3)] <- 1
}
TL4_1 <- cbind((matrix(0,18,dim(TL3_2)[2])),L4_1)
TL4_1.t <- cbind(TL4_1, (matrix(0,dim(TL4_1)[1],(897-dim(TL4_1)[2]) )))

#2000-2003
L4_2 <- diag(18*4)
TL4_2 <- cbind((matrix(0,dim(L4_2)[1],dim(TL4_1)[2])),L4_2)
TL4_2.t <- cbind(TL4_2, matrix(0,dim(TL4_2)[1],(897-dim(TL4_2)[2]) ))

#2004-2008
#2004
L5_1 <- matrix(0,25,18)
for(i in 1:13){
  L5_1[i,i] <- 1
}
L5_1[15,14] <- L5_1[17,15] <- L5_1[23,16] <- L5_1[24,17] <- L5_1[25,18] <- 1
TL5_1 <- cbind((matrix(0,25,dim(TL4_2)[2])),L5_1)
TL5_1.t <- cbind(TL5_1, (matrix(0,dim(TL5_1)[1],(897-dim(TL5_1)[2]) )))

#2005-2008
L5_2 <- diag(25*4)
TL5_2 <- cbind((matrix(0,dim(L5_2)[1],dim(TL5_1)[2])), L5_2)
TL5_2.t <- cbind(TL5_2, (matrix(0,dim(TL5_2)[1],(897-dim(TL5_2)[2]) )))

#2009-2021
#2009
L6_1 <- matrix(0,27,25)
for(i in 1:15){
  L6_1[i,i] <- 1
}
for(i in 16:25){
  L6_1[(i+2),i] <- 1
}
TL6_1 <- cbind((matrix(0,27, dim(TL5_2)[2])),L6_1)
TL6_1.t <- cbind(TL6_1, (matrix(0, dim(TL6_1)[1], (897-dim(TL6_1)[2]))))

L6_2 <- diag(27*12) #from 2010 to 2021
TL6_2 <- cbind((matrix(0,27*12,dim(TL6_1)[2])),L6_2)
TL6_2.t <- cbind(TL6_2, (matrix(0, dim(TL6_2)[1],(897-dim(TL6_2)[2]))))


TL <- rbind((rbind((rbind((rbind((rbind((rbind((rbind((rbind((rbind((rbind(TL1.t,TL2_1.t)),TL2_2.t)),TL3_1.t)),TL3_2.t)),TL4_1.t)),TL4_2.t)),TL5_1.t)),TL5_2.t)),TL6_1.t)),TL6_2.t)
dim(TL)
#big <- as.data.frame(TL)

#
capdistinv.wide <- as.matrix(mematrix1975_2021)
capdistinv.wide.resc <- capdistinv.wide*100 #re-scaled

sacL.resc$rho

#Identity matrix
I <- diag(897) 

Z1 <- (I - coef(sacL.resc)[4]*TL - coef(sacL.resc)[1]*capdistinv.wide.resc) 
I_Z1 <- solve(Z1)

#SR
indirect_SR_1975 = mean(rowSums(I_Z1[871:897,871:897])-diag(I_Z1[871:897,871:897])) 
indirect_SR_1975 #SR in Model 1 (Table 3 in Suppl materials)
#[1] -0.01563868

#or
I <- diag(27) # create identity matrix (N = 27)
matrix27x27 = capdistinv.wide.resc[871:897,871:897]
fx <- solve(I - (coef(sacL.resc)[1]*matrix27x27))
indirect_SR_1975 <- mean(rowSums(fx)-diag(fx))
indirect_SR_1975
#[1] -0.01563868

#LRSS = (I - phi * I - rho * W), based on Cook et al. (2023)
I <- diag(27) # N = 27
matrix27x27 = capdistinv.wide.resc[871:897,871:897]

matrix = solve(I - coef(sacL.resc)[4] * I - coef(sacL.resc)[1] * matrix27x27)
lrss1975 = mean(rowSums(matrix)-diag(matrix))  
lrss1975 #LRSS in Model 1 (Table 3 in Suppl materials) 
#[1] -0.504933

####################################################################################

##SR and LRSS effects (1991-2021) for Model 2 (Table 3)
####################################################################################
#
#building L multiplier (for unbalanced data)
#1991-2021

#1991-1996: 14 Nato countries (we don't have gdp data for Can until 1997)
#1997-1998: 15 (with Can)
#1999-2003: 18
dim(data.ordered1991)

#
L1 <- diag(14*5) #we leave the first year as containing only 0
TL1 <- rbind((matrix(0,14,dim(L1)[2])),L1)
TL1.t <- cbind(TL1,(matrix(0,dim(TL1)[1],(680 -dim(TL1)[2])))) #680 because our dataset has as many col.s

#
L2_1 <- matrix(0,15,14) #adding Can for 1997-1998
L2_1[1,1] <- 1 #US

for(i in 3:15){
  L2_1[i,(i-1)] <- 1
}
TL2_1 <- cbind((matrix(0,15,dim(TL1)[2])),L2_1)
TL2_1.t <- cbind(TL2_1,(matrix(0,dim(TL2_1)[1],(680-dim(TL2_1)[2]) )))

L2_2 <- diag(15) #year 1998
TL2_2 <- cbind((matrix(0,15,dim(TL2_1)[2])),L2_2)
TL2_2.t <- cbind(TL2_2,(matrix(0,dim(TL2_2)[1],(680-dim(TL2_2)[2]) )))

#1999-2003
#1999
L3_1 <- matrix(0,18,15) #18 - with 3 additional members, 15 - from the previous year
for(i in 1:10){
  L3_1[i,i] <- 1
}
for(i in 14:18){
  L3_1[i,(i-3)] <- 1
}
TL3_1 <- cbind((matrix(0,18,dim(TL2_2)[2])),L3_1)
TL3_1.t <- cbind(TL3_1, (matrix(0,dim(TL3_1)[1],(680-dim(TL3_1)[2]) )))

#2000-2003
L3_2 <- diag(18*4)
TL3_2 <- cbind((matrix(0,dim(L3_2)[1],dim(TL3_1)[2])),L3_2)
TL3_2.t <- cbind(TL3_2, matrix(0,dim(TL3_2)[1],(680-dim(TL3_2)[2]) ))

#2004-2008
#2004
L4_1 <- matrix(0,25,18)
for(i in 1:13){
  L4_1[i,i] <- 1
}
L4_1[15,14] <- L4_1[17,15] <- L4_1[23,16] <- L4_1[24,17] <- L4_1[25,18] <- 1
TL4_1 <- cbind((matrix(0,25,dim(TL3_2)[2])),L4_1)
TL4_1.t <- cbind(TL4_1, (matrix(0,dim(TL4_1)[1],(680-dim(TL4_1)[2]) )))

#2005-2008
L4_2 <- diag(25*4)
TL4_2 <- cbind((matrix(0,dim(L4_2)[1],dim(TL4_1)[2])), L4_2)
TL4_2.t <- cbind(TL4_2, (matrix(0,dim(TL4_2)[1],(680-dim(TL4_2)[2]) )))

#2009-2021
#2009
L5_1 <- matrix(0,27,25)
for(i in 1:15){
  L5_1[i,i] <- 1
}
for(i in 16:25){
  L5_1[(i+2),i] <- 1
}
TL5_1 <- cbind((matrix(0,27, dim(TL4_2)[2])),L5_1)
TL5_1.t <- cbind(TL5_1, (matrix(0, dim(TL5_1)[1], (680-dim(TL5_1)[2]))))

L5_2 <- diag(27*12) #from 2010 to 2021
TL5_2 <- cbind((matrix(0,27*12,dim(TL5_1)[2])),L5_2)
TL5_2.t <- cbind(TL5_2, (matrix(0, dim(TL5_2)[1],(680-dim(TL5_2)[2]))))

TL <- rbind((rbind((rbind((rbind((rbind((rbind((rbind((rbind(TL1.t,TL2_1.t)),TL2_2.t)),TL3_1.t)),TL3_2.t)),TL4_1.t)),TL4_2.t)),TL5_1.t)),TL5_2.t)
dim(TL)
#Identity matrix
I <- diag(680) 

dim(capdistinv.wide1991.resc) #from previously

Z1 <- (I - coef(sacL.resc1991)[4]*TL - coef(sacL.resc1991)[1]*capdistinv.wide1991.resc) 
I_Z1 <- solve(Z1)

#SR
indirect_SR_1991 = mean(rowSums(I_Z1[654:680,654:680])-diag(I_Z1[654:680,654:680]))  
indirect_SR_1991 #SR in Model 2 (Table 3 in Suppl materials) 
#[1] -0.02008391 

#or
I <- diag(27) # create identity matrix (N = 27)
matrix27x27.resc <- capdistinv.wide1991.resc[654:680,654:680]

fx <- solve(I - (coef(sacL.resc1991)[1]*matrix27x27.resc))
indirect_SR_1991 <- mean(rowSums(fx)-diag(fx))
indirect_SR_1991
#[1] -0.02008391

#
#LRSS = (I - phi * I - rho * W)
I <- diag(27) # create identity matrix (N = 27)
matrix27x27 = capdistinv.wide1991.resc[654:680,654:680]

matrix = solve(I - coef(sacL.resc1991)[4] * I - coef(sacL.resc1991)[1] * matrix27x27)
lrss1991 = mean(rowSums(matrix)-diag(matrix))  
lrss1991 #LRSS in Model 2 (Table 3 in Suppl materials)
#[1] -0.4198715

####################################################################################

##SR and LRSS effects (1999-2021) for Model 3 (Table 3)
####################################################################################

#building L multiplier (for unbalanced data)
#18 members - from 1999 to 2003
L1 <- diag(18*4)
TL1 <- rbind((matrix(0,18,dim(L1)[2])),L1)
TL1.t <- cbind(TL1,(matrix(0,dim(TL1)[1],(566 - dim(TL1)[2]))))

#25 members - from 2004 until 2008
L2_1 <- matrix(0,25,18)
for(i in 1:13){
  L2_1[i,i] <- 1
}
L2_1[15,14] <- L2_1[17,15] <- L2_1[23,16] <- L2_1[24,17] <- L2_1[25,18] <- 1
TL2_1 <- cbind((matrix(0,25,dim(TL1)[2])),L2_1) 
TL2_1.t <- cbind(TL2_1, (matrix(0,dim(TL2_1)[1],(566 - dim(TL2_1)[2]))))

L2_2 <- diag(25*4) #from 2005 to 2008
TL2_2 <- cbind((matrix(0,25*4,dim(TL2_1)[2])),L2_2)
TL2_2.t <- cbind(TL2_2,(matrix(0,dim(TL2_2)[1],(566 - dim(TL2_2)[2]))))

#27 members from 2009 to 2021
L3_1 <- matrix(0,27,25)
for(i in 1:15){
  L3_1[i,i] <- 1
}
for(i in 16:25){
  L3_1[(i+2),i] <- 1
}
TL3_1 <- cbind((matrix(0,27, dim(TL2_2)[2])),L3_1)
TL3_1.t <- cbind(TL3_1, (matrix(0, dim(TL3_1)[1], (566-dim(TL3_1)[2]))))

L3_2 <- diag(27*12) #from 2010 to 2021
TL3_2 <- cbind((matrix(0,27*12,dim(TL3_1)[2])),L3_2)
TL3_2.t <- cbind(TL3_2, (matrix(0, dim(TL3_2)[1],(566-dim(TL3_2)[2]))))

TL <- rbind((rbind((rbind((rbind(TL1.t,TL2_1.t)),TL2_2.t)),TL3_1.t)),TL3_2.t)
dim(TL)
#Identity matrix
I <- diag(566) 

Z1 <- (I - coef(sacL.resc1999)[4]*TL - coef(sacL.resc1999)[1]*capdistinv.wide1999.resc) 
I_Z1 <- solve(Z1)

#SR
indirect_SR_1999 = mean(rowSums(I_Z1[540:566,540:566])-diag(I_Z1[540:566,540:566])) 
indirect_SR_1999 #SR in Model 3 (Table 3 in Suppl materials)
#[1] -0.02921758

#or
I <- diag(27) # create identity matrix (N = 27)
matrix27x27 = capdistinv.wide1999.resc[540:566,540:566]
fx <- solve(I - (coef(sacL.resc1999)[1]*matrix27x27))
indirect_SR_1999 <- mean(rowSums(fx)-diag(fx))
indirect_SR_1999
#[1] -0.02921758

#
#LRSS = (I - phi * I - rho * W)
I <- diag(27) # create identity matrix (N = 27)
matrix27x27 = capdistinv.wide1999.resc[540:566,540:566]

matrix = solve(I - coef(sacL.resc1999)[4] * I - coef(sacL.resc1999)[1] * matrix27x27)
lrss1999 = mean(rowSums(matrix)-diag(matrix))  
lrss1999 #LRSS in Model 3 (Table 3 in Suppl materials)
#[1] -0.4770264

list.indirect_SR = c(indirect_SR_1975, indirect_SR_1991, indirect_SR_1999)
list.lrss = c(lrss1975, lrss1991, lrss1999)

m2 = matrix(NA, 2, 4)
colnames(m2) <- c("", "1975-2021", "1991-2021", "1999-2021")
m2[,1] = c("indirect", "LRSS")
m2[1,2:4] = list.indirect_SR
m2[2,2:4] = list.lrss
m2.df=as.data.frame(m2)

write_xlsx(m2.df, "Table3_Spatial effects_SM.xlsx")

####################################################################################

